local<-1
setwd("./output")
  
typename<-"stationary"


obsdraws<- as.data.table(read.dta13(paste0("psid_obs",typename,".dta")))
obsdraws_cons<-as.data.table(read.dta13(paste0("psid_obs_cons",typename,".dta")))
#NEED TO CORRECT DATAMOMENTS AND INVVARCOV:
temp<-read.table(paste0("wagemomentsstationary_jump.txt"),row.names=1,fill=TRUE)
temp<-temp[rownames(temp)%in%c("prodtype_1","prodtype_2","prodtype_3","age","age2","mod","sev","old","mod_old","sev_old"),]
wageregmoments<-temp[,1]
wageregse<-as.numeric(temp[,2])
temp<-read.table(paste0("spwagemomentsstationary_jump.txt"),row.names=1,fill=TRUE)
temp<-temp[rownames(temp)%in%c("prodtype_1","prodtype_2","prodtype_3","sp_age","sp_age2","sp_mod","sp_old","sp_mod_old"),]
spwageregmoments<-temp[,1]
spwageregse<-as.numeric(temp[,2])

temp<-read.table(paste0("wagevarsstationary_jump.txt"),row.names=1,fill=TRUE)
wagevarmoments<-temp[,1]
wagevarse<-temp[,2]


#invvarcov<-solve(as.matrix(read.table("~/Dropbox/DynamicDI/Low and Pistaferri/Code LP AER_Simulations/Code HL/varcovar30.inp",row.names=NULL)))
temp<-as.data.table(read.table("consmoments_all.txt"))
temp$V2<-as.numeric(as.character(temp$V2))
temp$V3<-as.numeric(as.character(temp$V3))
#dropping PT work for spouses:
temp<-temp[!grepl("sp_part_part",temp$V1),]
temp<-temp[!grepl("part_part_",temp$V1)]
consmoments_all<-temp[,V2] 
consse_all<-temp[,V3] 
temp<-as.data.table(read.table("consmoments_single.txt"))
temp<-temp[!grepl("part_part_",temp$V1)]
temp$V2<-as.numeric(as.character(temp$V2))
temp$V3<-as.numeric(as.character(temp$V3))
consmoments_single<-temp[,V2]
consse_single<-temp[,V3]  
temp<-as.data.table(read.table("consmoments_married.txt"))
#dropping PT work for spouses:
temp<-temp[!grepl("sp_part_part",temp$V1),]
temp$V2<-as.numeric(as.character(temp$V2))
temp$V3<-as.numeric(as.character(temp$V3))
consmoments_married<-temp[,V2] 
consse_married<-temp[,V3] 


temp<-as.data.table(read.table("consmoments_simple_all.txt"))
temp$V2<-as.numeric(as.character(temp$V2))
temp$V3<-as.numeric(as.character(temp$V3))
#dropping PT work for spouses:
temp<-temp[!grepl("sp_part_part",temp$V1),]
temp<-temp[!grepl("part_part",temp$V1)]
consmoments_simple_all<-temp[,V2] 
consse_simple_all<-temp[,V3] 

temp<-as.data.table(read.table("consmoments_simple_all.txt"))
temp<-temp[!grepl("part_part_",temp$V1)]
temp<-temp[!grepl("sp_",temp$V1)]
temp$V2<-as.numeric(as.character(temp$V2))
temp$V3<-as.numeric(as.character(temp$V3))
consmoments_simple_single<-temp[,V2]
consse_simple_single<-temp[,V3]  
temp<-as.data.table(read.table("consmoments_simple_all.txt"))
#dropping PT work for spouses:
temp<-temp[!grepl("sp_part_part",temp$V1),]
temp$V2<-as.numeric(as.character(temp$V2))
temp$V3<-as.numeric(as.character(temp$V3))
consmoments_simple_married<-temp[,V2] 
consse_simple_married<-temp[,V3] 


temp<-read.table("workmoments_agg.txt")
temp<-temp[!grepl("E=1/",rownames(temp)),]
workmoments_simple<-temp[,1]
workse_simple<-temp[,2] 

temp<-read.table("workmoments_1.txt")
workfullmoments_married<-temp[!grepl("E=1/",rownames(temp)),1]
workfullse_married<-temp[!grepl("E=1/",rownames(temp)),2]
workpartmoments_married<-temp[grepl("E=1/",rownames(temp)),1]
workpartse_married<-temp[grepl("E=1/",rownames(temp)),2]

temp<-temp[!grepl("E=1/",rownames(temp)),]
workmoments_mar<-temp[,1]
workse_mar<-temp[,2] 
#Dropping PT work:
#workmoments_mar<-temp[1:6,1]
#workse_mar<-temp[1:6,2] 
temp<-read.table("workmoments_0.txt")
workfullmoments_single<-temp[!grepl("E=1/",rownames(temp)),1]
workfullse_single<-temp[!grepl("E=1/",rownames(temp)),2]
workpartmoments_single<-temp[grepl("E=1/",rownames(temp)),1]
workpartse_single<-temp[grepl("E=1/",rownames(temp)),2]

temp<-temp[!grepl("E=1/",rownames(temp)),]
workmoments_sing<-temp[,1]
workse_sing<-temp[,2] 
# workmoments_sing<-temp[1:6,1]
# workse_sing<-temp[1:6,2] 
temp<-read.table("spworkmoments_1.txt")
spworkfullmoments<-temp[!grepl("E=1/",rownames(temp)),1]
spworkfullse<-temp[!grepl("E=1/",rownames(temp)),2]
spworkpartmoments<-temp[grepl("E=1/",rownames(temp)),1]
spworkpartse<-temp[grepl("E=1/",rownames(temp)),2]

temp<-read.table("workmoments_agg.txt")
temp<-temp[!grepl("E=1/",rownames(temp)),]
workfullmoments_simple<-temp[,1]
workfullse_simple<-temp[,2] 


temp<-read.table("stockDImoments_agg.txt")
stockDImoments_simple<-temp$EST[]
stockDIse_simple <- temp$SE[]

temp<-read.table("stockDImoments_work.txt")
temp<-read.table("stockDImoments.txt")
stockDImoments<-temp$EST[]
DIse <- temp$SE[]
stockDImoments_all<-temp$EST[]
stockDIse_all <- temp$SE[]
stockDImoments_single<-temp$EST[grepl("single",rownames(temp))]
stockDIse_single <- temp$SE[grepl("single",rownames(temp))]
stockDImoments_married<-temp$EST[!grepl("single",rownames(temp))]
stockDIse_married <- temp$SE[!grepl("single",rownames(temp))]

temp<-read.table("compDImoments_agg.txt")
compDImoments_simple<-temp$EST[]
compDIse_simple<-temp$SE[]

#stockDImoments<-temp$EST[-c(1,4,7,10)]
#DIse <- temp$SE[-c(1,4,7,10)]
temp<-read.table("compDImoments.txt")
compDImoments<-temp$EST[]
compDIse<-temp$SE[]
compDImoments_all<-temp$EST[]
compDIse_all <- temp$SE[]
compDImoments_single<-temp$EST[grepl("single",rownames(temp))]
compDIse_single <- temp$SE[grepl("single",rownames(temp))]
compDImoments_married<-temp$EST[!grepl("single",rownames(temp))]
compDIse_married <- temp$SE[!grepl("single",rownames(temp))]

#compDImoments<-temp$EST[-c(5,6,11,12)]
#compDIse<-temp$SE[-c(5,6,11,12)]

setwd("../eventstudy/PSID_eventstudy")

rform_events_mod<-as.data.table(read.csv("./mod_eventmeans_table.csv"))
rform_events_sev<-as.data.table(read.csv("./sev_eventmeans_table.csv"))
eventDImoments_single<-c(unlist(rform_events_mod[X%in%c("as.factor(married)0"),'Estimate']),
                         unlist(rform_events_sev[X%in%c("as.factor(married)0"),'Estimate']))
eventDIse_single<-c(unlist(rform_events_mod[X%in%c("as.factor(married)0"),'Std..Error']),
                    unlist(rform_events_sev[X%in%c("as.factor(married)0"),'Std..Error']))
eventDImoments_married<-c(unlist(rform_events_mod[X%in%c("as.factor(married)1"),'Estimate']),
                          unlist(rform_events_sev[X%in%c("as.factor(married)1"),'Estimate']))
eventDIse_married<-c(unlist(rform_events_mod[X%in%c("as.factor(married)1"),'Std..Error']),
                     unlist(rform_events_sev[X%in%c("as.factor(married)1"),'Std..Error']))
#single then married:
eventDImoments<-c(eventDImoments_single,eventDImoments_married)
eventDIse<-c(eventDIse_single,eventDIse_married)

rform_events_mod<-as.data.table(read.csv("./mod_agg_eventmeans_table.csv"))
rform_events_sev<-as.data.table(read.csv("./sev_agg_eventmeans_table.csv"))
eventDImoments_simple<-c(unlist(rform_events_mod[X%in%c("Estimate"),'x']),
                         unlist(rform_events_sev[X%in%c("Estimate"),'x']))
eventDIse_simple<-c(unlist(rform_events_mod[X%in%c("Std. Error"),'x']),
                    unlist(rform_events_sev[X%in%c("Std. Error"),'x']))

setwd("../../output")

temp<-read.table("flowDImoments_agg.txt")
flowDImoments_simple<-temp$EST[]
flowDIse_simple<- temp$SE[]


temp<-read.table("flowDImomentsbymar.txt")
flowDImoments<-temp$EST
names(flowDImoments)<-rownames(temp)
flowDIse<- temp$SE
names(flowDIse)<-rownames(temp)
flowDImoments_single<-temp$EST[grepl("M=0",rownames(temp))]
flowDIse_single <- temp$SE[grepl("M=0",rownames(temp))]
flowDImoments_married<-temp$EST[!grepl("M=0",rownames(temp))]
flowDIse_married <- temp$SE[!grepl("M=0",rownames(temp))]

temp<-read.table("contflowDImomentsbymar.txt")
contflowDImoments<-temp$EST
names(contflowDImoments)<-rownames(temp)
contflowDIse<- temp$SE
names(contflowDIse)<-rownames(temp)
contflowDImoments_single<-temp$EST[grepl("M=0",rownames(temp))]
contflowDIse_single <- temp$SE[grepl("M=0",rownames(temp))]
contflowDImoments_married<-temp$EST[!grepl("M=0",rownames(temp))]
contflowDIse_married <- temp$SE[!grepl("M=0",rownames(temp))]


temp<-read.table("jointworkmoments.txt")
jointworkmoments<-temp$EST #(no work, no spwork), (no work, sp work), (work, no spwork), (work, spwork)
names(jointworkmoments)<-rownames(temp)
jointworkse<-temp$SE
names(jointworkse)<-rownames(temp)
temp<-read.table("jointworkpartmoments.txt")
jointworkpartmoments<-temp$EST #(no work, no spwork), (no work, sp work), (work, no spwork), (work, spwork)
names(jointworkpartmoments)<-rownames(temp)
jointworkpartse<-temp$SE
names(jointworkpartse)<-rownames(temp)

#Loading wage parameters from a regression without selection correction (for post-estimation evaluation of fit):
wagevals_noselection<-read.table(paste0("wageparams_noselection",typename,".txt"),row.names=1,fill=TRUE)
spwagevals_noselection<-read.table(paste0("spwageparams_noselection",typename,".txt"),row.names=1,fill=TRUE)
wagevals_noselection<-wagevals_noselection[c("diffsev","diffmod","diffage","diffage2"),]
spwagevals_noselection<-spwagevals_noselection[c("diffsp_dis","diffsp_age","diffsp_age2"),]
wagenames<-c("$H=2$","$H=1$","$Age$","$Age^2/100$","Married")
spwagenames<-c("$H=1$","$Age$","$Age^2/100$")

#NOTE: need to adjust more than just this parameter if changing number of house types
numtypes<-3
agelength<-1

particles<-5 # no. particles in the particle swarm search
partic_mult<-1
#print(mpi.universe.size()-1)
seed=6327


numS=10
numsims<-read.table(paste0("numtypes",typename,".txt"))
numsims<-numsims[!rownames(numsims)=="Total",]*numS


testmar<-function(time){
  out<-NULL
  for(ty in 1:length(marriageprob_init)){
    out[ty]<-marriageprob_init[ty]
    for(t in 1:time){
      out[ty]<-out[ty]*marriageprob[[ty]][t,1]+(1-out[ty])*marriageprob[[ty]][t,2]
    }
  }
  return(out)
}

testhealth<-function(time){
  out<-matrix(c(1,0,0,0,0,0),3,6,byrow=TRUE)
  for(ty in 1:3){
    for(t in 1:time){
      for(s in 1:6){
        out[ty,s]<-out[ty,]%*%healthprob[[ty]][t,,s]
      }
    }
  }
  return(out)
}






if(modeltype=="simplified"){
  health_simpleprefs=TRUE
  single_sameprefs=TRUE
  single_sameallowanceprob=TRUE
  moments_simple=TRUE
} else{
  health_simpleprefs=FALSE
  single_sameprefs=FALSE
  single_sameallowanceprob=FALSE
  moments_simple=FALSE
}
ptwork_halfcost=FALSE
single_samecost=FALSE
ptwork_halfdis=FALSE
single_samedis=FALSE
spouse_samejobrisk=FALSE
work_noallowanceprob=FALSE


fullpar<-NULL
for(particle in 1:(partic_mult*particles)){
  if(local==0) guessvals<-read.csv(paste0('../../particlefiles/outputvals_fixmar',typename,particle,'.csv'),header=FALSE)
  if(local==1) guessvals<-read.csv(paste0('../../StructuralModel/particlefiles/outputvals_fixmar',typename,particle,'.csv'),header=FALSE)
  guessvals<-apply(guessvals,MARGIN=2, function(x) as.numeric(as.character(x)))
  guessvals<-guessvals[,-1]
  guessvals<-guessvals[guessvals[,ncol(guessvals)]==999,]
  guessvals<-guessvals[,-ncol(guessvals)]
  newguess<-unlist(as.vector(guessvals[which.min(guessvals[,dim(guessvals)[2]]),]))
  fullpar<-cbind(fullpar,newguess)
}
newguess<-fullpar[-dim(fullpar)[1],which.min(fullpar[dim(fullpar)[1],])]
fullpar<-fullpar[-dim(fullpar)[1],]
names(newguess)<-c("mod","sev","age","age2","Type1","Type2","Type3",
                   "spmod","spage","spage2","spType1","spType2","spType3",
                   "wagevar","spwagevar","wagecorr","noisevar","spnoisevar",
                   "theta_mod","theta_sev",
                   "spmu","sptheta",
                   "eta","etapart","eta_single","etapart_single","eta_mod","eta_sev",
                   "speta","spetapart","speta_mod",
                   "delta","delta_single", "spdelta",
                   "F1","F2","F3","Fpart1","Fpart2","Fpart3","Fold","Fpartold","F1_single","F2_single","F3_single","Fpart1_single","Fpart2_single",
                   "Fpart3_single","Fold_single","Fpartold_single","spF1","spF2","spFpart1","spFpart2","spFold","spFpartold",
                   "pi0.young.married","pi1.young.married","pi2.young.married","pi0.old.married","pi1.old.married","pi2.old.married",
                   "pi0.young.single","pi1.young.single","pi2.young.single","pi0.old.single","pi1.old.single","pi2.old.single",
                   "pi0.young.married.work","pi1.young.married.work","pi2.young.married.work","pi0.old.married.work","pi1.old.married.work","pi2.old.married.work",
                   "pi0.young.single.work","pi1.young.single.work","pi2.young.single.work","pi0.old.single.work","pi1.old.single.work","pi2.old.single.work")

upper<-c(0.05,0.05,0.2,0.2,3,3,5, #Ref person work params
         0.05,0.2,0.2,1,1,2, #spouse work params
         c(1,1,0.5,1,1), #work variance params
         0,0, #Thetas: disutility of health
         1.5,0.2, #spmu and sptheta: disutility of marriage and spousal health
         0,0, #etas: disutility of work, married
         0,0, #etas: disutility of work, single
         0, 0, #etas: disutility of work, when disabled
         0.2,0,0, #etas: disutility of work, for spouse
         0.2,0.6,1, #deltas: job loss probabilities
         0,1,3, #F: fixed cost of full-time work, married
         0,0,0, #F: fixed cost of part-time work, married
         #         0.5,2,4, #F: fixed cost of part-time work, married
         0.5,0, #F: extra fixed cost of work for being old, married
         #0.5,2, #F: extra fixed cost of work for being old, married
         0,2,3, #F: fixed cost of full-time work, single
         0,0,0, #F: fixed cost of part-time work, single
         #         0.5,2,3, #F: fixed cost of part-time work, single
         0.5,0, #F: extra fixed cost of work for being old, single
         #         0.5,2, #F: extra fixed cost of work for being old, single
         0,2, #F: fixed costs of full-timework for spouse
         0,0, #F: fixed costs of part-timework for spouse
         #         3,3, #F: fixed costs of part-timework for spouse
         0.5,0, #F: extra fxed cost of work for being old, spouse
         #0.5,1, #F: extra fxed cost of work for being old, spouse
         0.4,1,1,0.5,1,1, #allowance probs, married
         0.4,1,1,0.5,1,1, #allowance probs, single
         #         0.1,0.1,0.1,0.1,0.1,0.1, #allowance probs, married working
         #         0.1,0.1,0.1,0.1,0.1,0.1) #allowance probs, single working
         0,0,0,0,0,0, #allowance probs, married working
         0,0,0,0,0,0 #allowance probs, single working
) 

lower<-c(rep(-1,7), #Ref person work params
         rep(-1,6), #spouse work params
         rep(0.001,5), #work variance params
         -1.5,-1.5,
         0.2,-1,
         #         -0.5,-0.5,
         -0.5,0,
         #            -1,-0.5,
         -1,0,
         -0.2,-1,
         #-1,-1,-1,
         -1,0,-1,
         0,0,0,
         0,0,0,
         0,0,0,
         0,0,
         0,0,0,
         0,0,0,
         0,0,
         0,0,
         0,0,
         0,0,
         0,0,0.1,0,0,0.3,
         0,0,0,0,0,0,
         0,0,0,0,0,0,
         0,0,0,0,0,0
)
names(upper)<-names(newguess)
names(lower)<-names(newguess)

if(modeltype=="simplified"){#Killing much heterogeneity between married and single workers:
  lower[c("theta_mod","eta_mod", "spmu", "eta_single", "delta_single", "F2", "Fold", "F2_single", "F3_single", "Fold_single", "spFold",
          "pi0.young.single","pi1.young.single","pi2.young.single","pi0.old.single","pi1.old.single","pi2.old.single")]<-0
  upper[c("theta_mod","eta_mod", "spmu", "eta_single", "delta_single", "F2", "Fold", "F2_single", "F3_single", "Fold_single", "spFold",
          "pi0.young.single","pi1.young.single","pi2.young.single","pi0.old.single","pi1.old.single","pi2.old.single")]<-0
  
}

temp<-wagesetup_fit(newguess,0,0,wage_gridsize=10,spwage_gridsize=5#,filesuffix=typename
)
Wage<-temp$Wage
spWage<-temp$spWage
wagepars<-temp$wagepars
spwagepars<-temp$spwagepars
rm(temp)






temp<-paramsetup(newguess,wagepars,spwagepars,
                 ptwork_halfcost=ptwork_halfcost,
                 single_samecost=single_samecost,
                 ptwork_halfdis=ptwork_halfdis,
                 single_samedis=single_samedis,
                 spouse_samejobrisk=spouse_samejobrisk,
                 single_sameprefs = single_sameprefs,
                 health_simpleprefs = health_simpleprefs,
                 single_sameallowanceprob=single_sameallowanceprob,
                 work_noallowanceprob=work_noallowanceprob)
params<-temp$params
param_input<-temp$param_input
newguess<-newguess
rm(temp)


marriageprob_init<-as.data.table(read.csv(paste0("marriagemeans",typename,".csv"),row.names=NULL))
marriageprob_init<-marriageprob_init[,married]


#spwagepars$wagevar<-0
wagecov<-matrix(c(wagepars$wagevar^2,wagepars$wagecorr*(wagepars$wagevar*spwagepars$wagevar),
                  wagepars$wagecorr*(wagepars$wagevar*spwagepars$wagevar),spwagepars$wagevar^2),2,2)


states<-statesetup(dir="./",filesuffix=typename,
                   asset_gridsize=20,wage_gridsize=10)
Asset<-states$Asset
foodstamps<-states$foodstamps
numkids<-states$numkids
healthprob<-states$healthprob



marriageprob<-states$marriageprob
for(aa in 1:3){
  for(bb in 1:3){
    marriageprob[[aa]][[bb]][,1]<-1
    marriageprob[[aa]][[bb]][,2]<-0
  }
}
oecd<-states$oecd

set.seed(seed)
.Random.seed<-c(10407L,  -888780854L,   762756803L,   780773672L,  -352617207L,  1841665334L, -2063451585L)
#partnershocks<-exp(rmvnorm(50*sum(numsims),c(0,0),sigma=wagecov))
partnershocks<-rmvnorm(50*sum(numsims),c(0,0))
shocks.wageraw <- partnershocks[,1]
shocks.spwageraw <- partnershocks[,2]
shocks.wagenoise <- rnorm(50*sum(numsims),0,wagepars$noisevar)
shocks.spwagenoise <- rnorm(50*sum(numsims),0,spwagepars$noisevar)
shocks.health<- runif(50*sum(numsims))
shocks.allow<- runif(50*sum(numsims))
shocks.reassess<- runif(50*sum(numsims))
shocks.jobloss<- runif(50*sum(numsims))
shocks.spjobloss<- runif(50*sum(numsims))
shocks.marriage<- runif(50*sum(numsims))
shocks.wagereal <- exp(shocks.wageraw*wagepars$wagevar)
shocks.spwagereal <- exp((shocks.spwageraw*(1-wagepars$wagecorr^2)^0.5 + shocks.wageraw*wagepars$wagecorr)*spwagepars$wagevar)

goldtol <- params$goldsearch_tolerance
shocks.wage<-shocks.wagereal
shocks.spwage<-shocks.spwagereal

select<-drawselect(numsims)
selectC<-drawselectC(numsims)


noflow=FALSE
nocomp=FALSE
nostockDI=TRUE
noeventDI=FALSE
nocontflow=TRUE
onlyFT=1
workDI=0
noSSDI=0
nospmoments=FALSE
noworkDImoments=TRUE

setwd("../../")
